home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / lapack / dgeqpf.f < prev    next >
Text File  |  1996-07-19  |  7KB  |  221 lines
  1.       SUBROUTINE DGEQPF( M, N, A, LDA, JPVT, TAU, WORK, INFO )
  2. *
  3. *  -- LAPACK test routine (version 2.0) --
  4. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  5. *     Courant Institute, Argonne National Lab, and Rice University
  6. *     March 31, 1993
  7. *
  8. *     .. Scalar Arguments ..
  9.       INTEGER            INFO, LDA, M, N
  10. *     ..
  11. *     .. Array Arguments ..
  12.       INTEGER            JPVT( * )
  13.       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
  14. *     ..
  15. *
  16. *  Purpose
  17. *  =======
  18. *
  19. *  DGEQPF computes a QR factorization with column pivoting of a
  20. *  real M-by-N matrix A: A*P = Q*R.
  21. *
  22. *  Arguments
  23. *  =========
  24. *
  25. *  M       (input) INTEGER
  26. *          The number of rows of the matrix A. M >= 0.
  27. *
  28. *  N       (input) INTEGER
  29. *          The number of columns of the matrix A. N >= 0
  30. *
  31. *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
  32. *          On entry, the M-by-N matrix A.
  33. *          On exit, the upper triangle of the array contains the
  34. *          min(M,N)-by-N upper triangular matrix R; the elements
  35. *          below the diagonal, together with the array TAU,
  36. *          represent the orthogonal matrix Q as a product of
  37. *          min(m,n) elementary reflectors.
  38. *
  39. *  LDA     (input) INTEGER
  40. *          The leading dimension of the array A. LDA >= max(1,M).
  41. *
  42. *  JPVT    (input/output) INTEGER array, dimension (N)
  43. *          On entry, if JPVT(i) .ne. 0, the i-th column of A is permuted
  44. *          to the front of A*P (a leading column); if JPVT(i) = 0,
  45. *          the i-th column of A is a free column.
  46. *          On exit, if JPVT(i) = k, then the i-th column of A*P
  47. *          was the k-th column of A.
  48. *
  49. *  TAU     (output) DOUBLE PRECISION array, dimension (min(M,N))
  50. *          The scalar factors of the elementary reflectors.
  51. *
  52. *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)
  53. *
  54. *  INFO    (output) INTEGER
  55. *          = 0:  successful exit
  56. *          < 0:  if INFO = -i, the i-th argument had an illegal value
  57. *
  58. *  Further Details
  59. *  ===============
  60. *
  61. *  The matrix Q is represented as a product of elementary reflectors
  62. *
  63. *     Q = H(1) H(2) . . . H(n)
  64. *
  65. *  Each H(i) has the form
  66. *
  67. *     H = I - tau * v * v'
  68. *
  69. *  where tau is a real scalar, and v is a real vector with
  70. *  v(1:i-1) = 0 and v(i) = 1; v(i+1:m) is stored on exit in A(i+1:m,i).
  71. *
  72. *  The matrix P is represented in jpvt as follows: If
  73. *     jpvt(j) = i
  74. *  then the jth column of P is the ith canonical unit vector.
  75. *
  76. *  =====================================================================
  77. *
  78. *     .. Parameters ..
  79.       DOUBLE PRECISION   ZERO, ONE
  80.       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
  81. *     ..
  82. *     .. Local Scalars ..
  83.       INTEGER            I, ITEMP, J, MA, MN, PVT
  84.       DOUBLE PRECISION   AII, TEMP, TEMP2
  85. *     ..
  86. *     .. External Subroutines ..
  87.       EXTERNAL           DGEQR2, DLARF, DLARFG, DORM2R, DSWAP, XERBLA
  88. *     ..
  89. *     .. Intrinsic Functions ..
  90.       INTRINSIC          ABS, MAX, MIN, SQRT
  91. *     ..
  92. *     .. External Functions ..
  93.       INTEGER            IDAMAX
  94.       DOUBLE PRECISION   DNRM2
  95.       EXTERNAL           IDAMAX, DNRM2
  96. *     ..
  97. *     .. Executable Statements ..
  98. *
  99. *     Test the input arguments
  100. *
  101.       INFO = 0
  102.       IF( M.LT.0 ) THEN
  103.          INFO = -1
  104.       ELSE IF( N.LT.0 ) THEN
  105.          INFO = -2
  106.       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  107.          INFO = -4
  108.       END IF
  109.       IF( INFO.NE.0 ) THEN
  110.          CALL XERBLA( 'DGEQPF', -INFO )
  111.          RETURN
  112.       END IF
  113. *
  114.       MN = MIN( M, N )
  115. *
  116. *     Move initial columns up front
  117. *
  118.       ITEMP = 1
  119.       DO 10 I = 1, N
  120.          IF( JPVT( I ).NE.0 ) THEN
  121.             IF( I.NE.ITEMP ) THEN
  122.                CALL DSWAP( M, A( 1, I ), 1, A( 1, ITEMP ), 1 )
  123.                JPVT( I ) = JPVT( ITEMP )
  124.                JPVT( ITEMP ) = I
  125.             ELSE
  126.                JPVT( I ) = I
  127.             END IF
  128.             ITEMP = ITEMP + 1
  129.          ELSE
  130.             JPVT( I ) = I
  131.          END IF
  132.    10 CONTINUE
  133.       ITEMP = ITEMP - 1
  134. *
  135. *     Compute the QR factorization and update remaining columns
  136. *
  137.       IF( ITEMP.GT.0 ) THEN
  138.          MA = MIN( ITEMP, M )
  139.          CALL DGEQR2( M, MA, A, LDA, TAU, WORK, INFO )
  140.          IF( MA.LT.N ) THEN
  141.             CALL DORM2R( 'Left', 'Transpose', M, N-MA, MA, A, LDA, TAU,
  142.      $                   A( 1, MA+1 ), LDA, WORK, INFO )
  143.          END IF
  144.       END IF
  145. *
  146.       IF( ITEMP.LT.MN ) THEN
  147. *
  148. *        Initialize partial column norms. The first n elements of
  149. *        work store the exact column norms.
  150. *
  151.          DO 20 I = ITEMP + 1, N
  152.             WORK( I ) = DNRM2( M-ITEMP, A( ITEMP+1, I ), 1 )
  153.             WORK( N+I ) = WORK( I )
  154.    20    CONTINUE
  155. *
  156. *        Compute factorization
  157. *
  158.          DO 40 I = ITEMP + 1, MN
  159. *
  160. *           Determine ith pivot column and swap if necessary
  161. *
  162.             PVT = ( I-1 ) + IDAMAX( N-I+1, WORK( I ), 1 )
  163. *
  164.             IF( PVT.NE.I ) THEN
  165.                CALL DSWAP( M, A( 1, PVT ), 1, A( 1, I ), 1 )
  166.                ITEMP = JPVT( PVT )
  167.                JPVT( PVT ) = JPVT( I )
  168.                JPVT( I ) = ITEMP
  169.                WORK( PVT ) = WORK( I )
  170.                WORK( N+PVT ) = WORK( N+I )
  171.             END IF
  172. *
  173. *           Generate elementary reflector H(i)
  174. *
  175.             IF( I.LT.M ) THEN
  176.                CALL DLARFG( M-I+1, A( I, I ), A( I+1, I ), 1, TAU( I ) )
  177.             ELSE
  178.                CALL DLARFG( 1, A( M, M ), A( M, M ), 1, TAU( M ) )
  179.             END IF
  180. *
  181.             IF( I.LT.N ) THEN
  182. *
  183. *              Apply H(i) to A(i:m,i+1:n) from the left
  184. *
  185.                AII = A( I, I )
  186.                A( I, I ) = ONE
  187.                CALL DLARF( 'LEFT', M-I+1, N-I, A( I, I ), 1, TAU( I ),
  188.      $                     A( I, I+1 ), LDA, WORK( 2*N+1 ) )
  189.                A( I, I ) = AII
  190.             END IF
  191. *
  192. *           Update partial column norms
  193. *
  194.             DO 30 J = I + 1, N
  195.                IF( WORK( J ).NE.ZERO ) THEN
  196.                   TEMP = ONE - ( ABS( A( I, J ) ) / WORK( J ) )**2
  197.                   TEMP = MAX( TEMP, ZERO )
  198.                   TEMP2 = ONE + 0.05D0*TEMP*
  199.      $                    ( WORK( J ) / WORK( N+J ) )**2
  200.                   IF( TEMP2.EQ.ONE ) THEN
  201.                      IF( M-I.GT.0 ) THEN
  202.                         WORK( J ) = DNRM2( M-I, A( I+1, J ), 1 )
  203.                         WORK( N+J ) = WORK( J )
  204.                      ELSE
  205.                         WORK( J ) = ZERO
  206.                         WORK( N+J ) = ZERO
  207.                      END IF
  208.                   ELSE
  209.                      WORK( J ) = WORK( J )*SQRT( TEMP )
  210.                   END IF
  211.                END IF
  212.    30       CONTINUE
  213. *
  214.    40    CONTINUE
  215.       END IF
  216.       RETURN
  217. *
  218. *     End of DGEQPF
  219. *
  220.       END
  221.